#ifndef BTL_MATHS_FLOATCOMPARISON
#define BTL_MATHS_FLOATCOMPARISON

// Implements ulp-based floating point comparison,
// based on the explanation at:
// http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
//
// note: Boost's Maths library includes some special functions
// to find the ULPs between two floats, or move to the next or previous
// float; but they do things I don't want, like raising exceptions
// if presented with infinities on NaNs

namespace btl
{
namespace maths
{

template <typename T> struct IEEEFloatingPointTraits;

template <typename T>
typename IEEEFloatingPointTraits<T>::bool_type
   almost_equal(T a, T b,
      typename IEEEFloatingPointTraits<T>::int_type max_ulps
      = IEEEFloatingPointTraits<T>::DefaultUlpTolerance);

template <typename T>
typename IEEEFloatingPointTraits<T>::int_type
   ulp_distance(T a, T b);

template <typename T>
typename IEEEFloatingPointTraits<T>::bool_type
   almost_zero(T x,
      T tolerance
      = std::numeric_limits<T>::epsilon());

template <typename T>
typename IEEEFloatingPointTraits<T>::bool_type is_nan(T x);

template <typename T>
typename IEEEFloatingPointTraits<T>::bool_type is_finite(T x);

} // namespace maths
} // namespace btl

namespace btl
{

using maths::almost_equal;
using maths::almost_zero;
using maths::is_nan;
using maths::is_finite;

} // namespace btl

// ====================================================================
// === Implementation

namespace btl
{
namespace maths
{

template <typename T>
inline bool is_finite_bits(const typename IEEEFloatingPointTraits<T>::uint_type& bits)
{
   return
      ((bits & IEEEFloatingPointTraits<T>::ExponentBits) != IEEEFloatingPointTraits<T>::ExponentBits);
}

template <typename T>
inline bool is_nan_bits(const typename IEEEFloatingPointTraits<T>::uint_type& bits)
{
   // NaN has the exponent bits set, and at least one mantissa bit set
   // (therefore, if you mask off the top bit, the result must be strictly greater than
   //  just the exponent bits set; if it's equal then it's just an infinity; if it's
   //  less, then it's a valid finite number)
   return
      ((bits & ~IEEEFloatingPointTraits<T>::TopBit) > IEEEFloatingPointTraits<T>::ExponentBits);
}

// IEEEFloatingPointTraits<T>::bool_type used for SFINAE
template <typename T>
inline typename IEEEFloatingPointTraits<T>::bool_type is_finite(T x)
{
   return is_finite_bits(*reinterpret_cast<typename IEEEFloatingPointTraits<T>::uint_type*>(&x));
}

// IEEEFloatingPointTraits<T>::bool_type used for SFINAE
template <typename T>
inline typename IEEEFloatingPointTraits<T>::bool_type is_nan(T x)
{
   return is_nan_bits(*reinterpret_cast<typename IEEEFloatingPointTraits<T>::uint_type*>(&x));
}

// IEEEFloatingPointTraits<T>::bool_type used for SFINAE
template <typename T>
inline typename IEEEFloatingPointTraits<T>::bool_type
   almost_equal(T a, T b, typename IEEEFloatingPointTraits<T>::int_type max_ulps)
{
   typedef typename IEEEFloatingPointTraits<T>::uint_type uint_type;
   typedef typename IEEEFloatingPointTraits<T>::int_type int_type;
   uint_type& aui = *reinterpret_cast<uint_type*>(&a);
   uint_type& bui = *reinterpret_cast<uint_type*>(&b);

   // Infinities aren't close to anything except themselves
   if (  !is_finite_bits<T>(aui) && is_finite_bits<T>(bui)
      || is_finite_bits<T>(aui) && !is_finite_bits<T>(bui))
      return false;

   // IEEE says NaNs are unequal to everything (even themselves)
   if (is_nan_bits<T>(aui) || is_nan_bits<T>(bui))
      return false;

   // transform from sign-magnitude to two's-complement
   if (aui < 0) aui = (IEEEFloatingPointTraits<T>::TopBit - aui);
   if (bui < 0) bui = (IEEEFloatingPointTraits<T>::TopBit - bui);

   int_type& ai = *reinterpret_cast<int_type*>(&aui);
   int_type& bi = *reinterpret_cast<int_type*>(&bui);

   int_type ulps = std::abs(bi - ai);
   return (ulps <= max_ulps);
}

// IEEEFloatingPointTraits<T>::int_type used for SFINAE
template <typename T>
typename IEEEFloatingPointTraits<T>::int_type ulp_distance(T a, T b)
{
   typedef typename IEEEFloatingPointTraits<T>::uint_type uint_type;
   typedef typename IEEEFloatingPointTraits<T>::int_type int_type;
   uint_type& aui = *reinterpret_cast<uint_type*>(&a);
   uint_type& bui = *reinterpret_cast<uint_type*>(&b);

   // transform from sign-magnitude to two's-complement
   if (aui < 0) aui = (IEEEFloatingPointTraits<T>::TopBit - aui);
   if (bui < 0) bui = (IEEEFloatingPointTraits<T>::TopBit - bui);

   int_type& ai = *reinterpret_cast<int_type*>(&aui);
   int_type& bi = *reinterpret_cast<int_type*>(&bui);

   return (bi - ai);
}

// IEEEFloatingPointTraits<T>::bool_type used for SFINAE
template <typename T>
typename IEEEFloatingPointTraits<T>::bool_type almost_zero(T x, T tolerance)
{
   return (std::abs(x) < tolerance);
}

template <typename T>
struct IEEEFloatingPointTraits {};

template <>
struct IEEEFloatingPointTraits<double>
{
   typedef bool bool_type;

   typedef long long int_type;
   typedef unsigned long long uint_type;

   static const uint_type TopBit
      = static_cast<uint_type>(1) << (sizeof(double)*8 - 1);
   static const uint_type ExponentBits
      = (~static_cast<uint_type>(0) << std::numeric_limits<double>::digits) & ~TopBit;
   static const uint_type MantissaBits
      = ~TopBit & ~ExponentBits;

   static const int_type DefaultUlpTolerance
      = 16;
};

template <>
struct IEEEFloatingPointTraits<float>
{
   typedef int ulp_int_type;
   typedef int int_type;
   typedef unsigned int uint_type;

   static const uint_type TopBit
      = uint_type(1) << (sizeof(float)*8 - 1);
   static const uint_type ExponentBits
      = (~uint_type(0) << std::numeric_limits<float>::digits) & ~TopBit;
   static const uint_type MantissaBits
      = ~TopBit & ~ExponentBits;

   static const int_type DefaultUlpTolerance
      = 4;
};

} // namespace maths
} // namespace btl

#endif // BTL_MATHS_FLOATCOMPARISON
